libname in "INSERT PATH TO SAS LIBRARY";
options linesize = 72;

data therms;
set in.therms;
array flip bush--reppty;
   do over flip;
   flip = 100 - flip;
   end;
run;

OPTIONS PAGESIZE = 58;
PROC IML WORKSIZE=300;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          This PROC IML routine calculates the line-of-            */
/*          sight measure of interobject similarity.  The            */
/*          reference for this procedure is:                         */
/*                                                                   */
/*                Rabinowitz, George (1976) A Procedure for          */
/*                Ordering Object Pairs Consistent with the          */
/*                Multidimensional Unfolding Model.                  */
/*                PSYCHOMETRIKA 41: 349-373.                         */
/*                                                                   */
/*          The proper citation for using this SAS procedure is:     */
/*                                                                   */
/*                Jacoby, William G. (1993) A SAS/IML Macro          */
/*                for Calculating the Line-of-Sight Measure          */
/*                of Interobject Dissimilarity.                      */
/*                PSYCHOMETRIKA 58: 511-512.                         */
/*                                                                   */
/*          This procedure assumes an input SAS data set that        */
/*          consists of N subjects' ratings of K objects.  Hence     */
/*          The input data are an N by K matrix.  There should       */
/*          be no missing values, and the scores should be           */
/*          set so that SMALLER values are more positive ratings,    */
/*          and HIGHER values are less positive ratings.             */
/*                                                                   */
/*-------------------------------------------------------------------*/

START NAME;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module NAME obtains the names of the variables           */
/*          from the input data set, by reading in a single          */
/*          observation, using the COLNAME specification, and        */
/*          then freeing the data.  The names of the rating          */
/*          variables are then stored in the row vector, NAMES.      */
/*          By default, the routine uses the most recently           */
/*          created SAS dataset for input data.  A different         */
/*          SAS dataset can be specified in the USE statement.       */
/*                                                                   */
/*-------------------------------------------------------------------*/

USE _LAST_;
READ POINT 1 VAR _NUM_ INTO DUMMY [COLNAME=NAMES];
FREE DUMMY;
FINISH;

START READIN;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module READIN reads the data from the input SAS          */
/*          dataset, and calculates the pairwise sums and            */
/*          differences from the rating scores.  Begin by            */
/*          looping through the columns of the input dataset.        */
/*                                                                   */
/*                I1:  Index for first var in current pair           */
/*             NAME1:  Name of first var in current pair             */
/*               ONE:  Vector of ratings of first object             */
/*                     in the current pair                           */       
/*                                                                   */
/*-------------------------------------------------------------------*/

   DO I1 = 1 TO NCOL(NAMES) - 1;
   NAME1 = NAMES[,I1];
   READ ALL VAR NAME1 INTO ONE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          I2:  Index for second var in current pair                */
/*       NAME2:  Name of second var in current pair                  */
/*                                                                   */
/*-------------------------------------------------------------------*/

     DO I2 = I1 + 1 TO NCOL(NAMES);
     NAME2 = NAMES[,I2];
     READ ALL VAR NAME2 INTO TWO;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Create vector of current sums (SUMVEC) and               */
/*          vector of current differences (DIFVEC).  Then            */
/*          sort SUMVEC from smallest to largest, and sort           */
/*          DIFVEC from largest to smallest.                         */
/*                                                                   */
/*-------------------------------------------------------------------*/

     SUMVEC = ONE + TWO;
     DIFVEC = ABS(ONE - TWO);
     FREE TWO;
     SORTVEC = SUMVEC;
     SUMVEC[RANK(SUMVEC),] = SORTVEC;
     SORTVEC = DIFVEC;
     DIFVEC[(NROW(DIFVEC) + 1) - RANK(DIFVEC),] = SORTVEC;
     FREE SORTVEC;
       IF I1 = 1 & I2 = 2 THEN DO;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          At first pair, create matrix of cumulative, sorted       */
/*          pairwise sums and differences.  NOBS is the total        */
/*          number of observations in the input data set, and        */
/*          NROWS is the maximum number of iterations that will      */
/*          be used in this run.  After the current sums and         */
/*          differences are concatenated to the SUM and DIF          */
/*          matrices, SUMVEC and DIFVEC can be freed.  At the        */
/*          end of the READIN module, both SUM and DIF will have     */
/*          NROWS rows and K(K-1)/2 columns.                         */
/*                                                                   */
/*-------------------------------------------------------------------*/

       NOBS = NROW(SUMVEC);
       NROWS = ROUND (NOBS / 25);
       SUM = CUSUM(SUMVEC[1:NROWS,]);
       DIF = CUSUM(DIFVEC[1:NROWS,]);
       FREE SUMVEC DIFVEC;
       END;
       IF I1 > 1 | I2 > 2 THEN DO;
       SUM = SUM || CUSUM(SUMVEC[1:NROWS,]);
       DIF = DIF || CUSUM(DIFVEC[1:NROWS,]);
       FREE SUMVEC DIFVEC;
       END;
     END;
   FREE ONE;
   END;
FINISH;

START CALC;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module CALC uses the cumulative, sorted SUM and          */
/*          DIF matrices to calculate the LOS matrix.                */
/*          Matrices are defined as follows:                         */
/*                                                                   */
/*          CURSTATS:  Vector of current LOS statistics              */
/*          STATNAME:  Vector of labels for printed output           */
/*            NOINCR:  Number of iterations needed to check          */
/*                     for no increase in ADEQU value                */
/*                                                                   */
/*-------------------------------------------------------------------*/       


   CURSTATS = J(1,3,0);
   STATNAME = {DISCRIM DENSE ADEQU};
   NOINCR = ROUND(NOBS / 150);
   IF NOINCR < 4 THEN NOINCR = 4;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Start looping through rows of pairwise sums and          */
/*          differences.  Sum corresponding entries together,        */
/*          obtain rank order of resultant sums.  Then use           */
/*          DISNUM and NDIST to obtain number of distinct            */
/*          ranks on this iteration.                                 */
/*                                                                   */
/*-------------------------------------------------------------------*/

   DO I = 1 TO NROWS;
   DISNUM = RANKTIE(SUM[I,] + DIF[I,]);
   NDIST = 0;
     DO J = 1 TO NCOL(DISNUM);
     IF DISNUM[1,J] > 0 THEN NDIST = NDIST + 1;
       IF J < NCOL(DISNUM) THEN DO;
         DO J2 = J+1 TO NCOL(DISNUM);
         IF DISNUM[1,J2] = DISNUM[1,J] THEN
            DISNUM[1,J2] = -1;
         END;
       END;
     END;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          CURSTATS contains the LOS statistics for the             */
/*          current row of the SUM and DIF matrices.  These          */
/*          values are concatenated onto the bottom of the           */
/*          STATS matrix, for later printing the iteration           */
/*          history.  CURMAX is the current maximum value of         */
/*          the ADEQU statistic.  The row in which this max          */
/*          ADEQU value occurs is indexed by the var CURLOC.         */
/*          The LOS statistics for the current max ADEQU value       */
/*          are contained in the vector MAXSTATS.  On each           */
/*          iteration, check the ADEQU for that row, and replace     */
/*          max ADEQU information, if necessary.                     */ 
/*                                                                   */
/*-------------------------------------------------------------------*/   

   CURSTATS[,1] = NDIST / (NCOL(DISNUM) - 1);
   CURSTATS[,2] = (NOBS - I) / (NOBS - 1);
   CURSTATS[,3] = CURSTATS[,1] # (CURSTATS[,2] ## 3);
     IF I = 1 THEN DO;
     STATS = CURSTATS;
     CURMAX = CURSTATS[,3];
     CURLOC = I;
     MAXSTATS = CURSTATS;
     END;
     IF I > 1 THEN STATS = STATS // CURSTATS;
     IF CURSTATS[,3] > CURMAX THEN DO;
     CURMAX = CURSTATS[,3];
     CURLOC = I;
     MAXSTATS=CURSTATS;
     END;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Check to see whether ADEQU has not increased over        */
/*          the last NOINCR rows.  If it has not increased, then     */
/*          terminate the procedure.                                 */
/*                                                                   */
/*-------------------------------------------------------------------*/

     IF I >= NOINCR THEN DO;
     C1 = STATS[(I-NOINCR+1):(I-1),3];
     C2 = STATS[(I-NOINCR+2):(I),3];
     IF SUM((C1 - C2) < 0) = 0 THEN GO TO TERM2;
     END;
   END;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          TERM1 is used if the procedure goes through all          */
/*          NROWS rows of the matrices.  The I value is              */
/*          decreased by 1, so it prints the correct number          */
/*          of iterations.                                           */
/*                                                                   */
/*-------------------------------------------------------------------*/

TERM1:
   I = I - 1;
   PRINT  'FINAL ITERATION NUMBER:' I,
     'ROUTINE TERMINATED BECAUSE IT REACHED MAX N OF ITERATIONS,,';
   GO TO DONE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          TERM2 is used if the ADEQU value converges.              */
/*                                                                   */
/*-------------------------------------------------------------------*/

TERM2:
   PRINT  'FINAL ITERATION NUMBER:' I,
     'ROUTINE TERMINATED BECAUSE ADEQUACY VALUE CONVERGED,,';
   GO TO DONE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*        DONE obtains the final ranking, based on combined sums     */
/*        and differences, and places it into the row vector,        */
/*        OBJRANK.  It also calculates the Spearman correlation      */
/*        between the ranking based on sums only, and that based     */
/*        on differences only (RHO).  It also prints out the         */
/*        iteration history and the LOS statistics for the           */
/*        final solution.                                            */
/*                                                                   */
/*-------------------------------------------------------------------*/

DONE:
   OBJRANK = RANKTIE(SUM[CURLOC,] + DIF[CURLOC,]);
   RHO = 1 - ((6 # SSQ(RANKTIE(SUM[CURLOC,]) - RANKTIE(DIF[CURLOC,])))
         /(NCOL(SUM) # ((NCOL(SUM)##2) - 1)));
   ITER = (1:NROW(STATS))`;   
   ITNAME = {'ITER #'};
   PRINT 'ITERATION HISTORY', ITER [COLNAME=ITNAME FORMAT=3.] 
     STATS [COLNAME = STATNAME FORMAT=7.3],,
     'FINAL LOS STATISTICS:',,
     'SELECTED OBJECT RANKING FROM ITERATION #:' CURLOC ,,
     MAXSTATS [COLNAME=STATNAME FORMAT=7.3] ,,
     'SUM-DIFF RHO:' RHO [FORMAT=7.3];
   FINISH;

START SHAPE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module SHAPE places the final LOS values, obtained       */
/*          from the OBJRANK vector, into a square                   */
/*          K by K matrix called DISSIM.                             */
/*                                                                   */
/*-------------------------------------------------------------------*/
                
DISSIM = J(NCOL(NAMES),NCOL(NAMES),0);
I = 0;
   DO I1 = 1 TO NCOL(NAMES) - 1;
     DO I2 = I1 + 1 TO NCOL(NAMES);
     I = I + 1;
     DISSIM[I1,I2] = OBJRANK[,I];
     DISSIM[I2,I1] = OBJRANK[,I];
     END;
   END;
PRINT 'FINAL MATRIX OF DISSIMILARITIES',,
      DISSIM [ROWNAME=NAMES COLNAME=NAMES];
FINISH;

START WRITE;

/*-------------------------------------------------------------------*/
/*                                                                   */
/*          Module WRITE creates a new SAS dataset, called           */
/*          DISSIM, from the previously-obtained DISSIM matrix.      */
/*          This dataset (with K observations and K variables)       */
/*          contains the pairwise LOS dissimilarity values.  The     */
/*          variable names are identical to those in the input       */
/*          SAS dataset.  Also, remember that the LOS values are     */
/*          dissimilarities-- lower values are more similar pairs    */
/*          and higher values are more dissimilar pairs              */
/*                                                                   */
/*-------------------------------------------------------------------*/

CREATE DISSIM FROM DISSIM [COLNAME = NAMES];
APPEND FROM DISSIM;
FINISH;

RUN NAME;
RUN READIN;
RUN CALC;
RUN SHAPE;
RUN WRITE;

PROC PRINT DATA=DISSIM;
RUN;
